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Abstract. The thermal structure of the post-shock region of a young supernova remnant 
(SNR) is heavily affected by two main physical effects, the back-reaction of accelerated cos- 
mic rays (CRs) and the Rayleigh-Taylor (RT) instabilities developing at the contact discon- 
tinuity between the ejecta and the shocked interstellar medium (ISM). Here, we investigate 
the role played by both physical mechanisms in the evolution of SNRs through detailed 3D 
MHD modeling. Our model describes the expansion of the remnant through a magnetized 
ISM, including consistently the initial ejecta clumping and the effects on shock dynamics 
due to back-reaction of accelerated CRs. We discuss the role of the initial ejecta dumpiness 
in developing strong instabilities at the contact discontinuity which may extend upstream to 
the main shock and beyond. 
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1. Introduction 

SNRs are the site where CR diffusive shock ac- 
celeration occurs. Indirect evidence of CRs ac- 
celeration comes from the separation between 
the forward shock and the contact disconti- 
nuity. In SN 1006, observations have shown 
that the azimuthal profile of the ratio between 
the forward shock and the contact discontinu- 
ity radii Rbw/Rcd is fairly uniform and much 
lower than predicted for a non-modified shock 
(Miceli et al. 2009). Recently Rakowski et al. 
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(2011) have found and analyzed clumps of 
ejecta close to or protruding beyond the main 
blast wave. These findings have been inter- 
preted as a consequence of the energy losses 
to CRs at the forward shock. However, sev- 
eral authors have shown that extreme energy 
losses to accelerate the CRs are needed to al- 
low a significant fraction of the ejecta to ap- 
proach or even overtake the forward shock (e.g 
Miceli et al. 2009). In addition, the observa- 
tions show that Rbw/Rcd is lower than predicted 
by non-modified shock models even in regions 
dominated by thermal emission where the CR 
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acceleration efficiency is supposed to be low 
(e.g. Miceli et al. 2009). 

An alternative in the data interpretation 
could be the ejecta structure of the explosion 
itself. In fact, density inhomogeneities in the 
ejecta can enhance the growth of RT instabili- 
ties potentially allowing them to reach the for- 
ward shock. The question is: can the ejecta 
clumping enhance the growth of RT instabil- 
ities up to a level that allows them to reach 
and possibly overtake the forward shock? In 
fact, there is a growing consensus that density 
clumping of ejecta is intrinsic at early phases 
of the remnant evolution. The density structure 
of ejecta therefore may naturally explain the 
small values of Rbw/Rcd observed in SNRs. 

Here we investigate this issue by develop- 
ing a 3D MHD model describing the expan- 
sion of a SNR through a magnetized medium, 
including consistently both the initial ejecta 
clumping and the effects on shock dynamics 
due to back-reaction of accelerated CRs. 

2. Three-dimensional MHD Modeling 

We improved the 3D MHD model discussed 
by Orlando et al. (2007, 2011). The SNR 
is modeled by numerically solving the time- 
dependent MHD equations of mass, momen- 
tum, and energy conservation, including con- 
sistently the initial ejecta clumping and the ef- 
fects on shock dynamics due to back-reaction 
of accelerated CRs (Orlando et al. 2011, in 
preparation). The calculations were performed 
using FLASH, an adaptive mesh refinement 
multiphysics code for astrophysical plasmas 
(Fryxell et al. 2000) extended with additional 
computational modules to handle the back- 
reaction of accelerated CRs. The latter process 
is included in the model by following the ap- 
proach of Ferrand et al. (2010) and extending 
their method to MHD models. In particular, our 
model includes an effective adiabatic index y e ff 
evolving in time: at each time-step of integra- 
tion the time- and space-dependent y e $ is cal- 
culated at the shock front and then is advected 
within the remnant, remaining constant in each 
fluid element. In addition to the approach of 
Ferrand et al. (2010) we allow y ei f to depend 
on local physical conditions and on the injec- 



tion rate of particles, by assuming that the in- 
dex depends on the obliquity angle between the 
upstream magnetic field and the normal to the 
shock. 



As for the density structure of the ejecta, 
we have investigated two different initial ejecta 
density profiles: the exponential profile that 
has been shown to be the most representative 
of explosion models for thermonuclear SNe 
(Dwarkadas & Chevalier 1998), and the power- 
law profile with n = 7 that has been used to 
represent deflagration models (Nomoto et al. 
1984). We assume that the initial ejecta have 
a clumpy structure, modeled by defining small 
clumps including a perturbation of the mass 
density. For each of the initial profiles, we ex- 
plore density clumps of ejecta with initial size 
either 2% or 4% of the initial remnant radius. 
Initially each clump has a random density con- 
trast derived from an exponential distribution 
that is characterized by a maximum density 
contrast. We explore maximum contrasts rang- 
ing between 1.5 and 5. 



As for the reference physical system, we 
adopted parameters appropriate to describe 
SN 1006: density of ISM n ism = 0.05 cnr 3 , en- 
ergy of the SN explosion £sn — 1.5 xlO 51 ergs, 
mass of ejecta M e j = 1 .4 M , shock velocity 
after 1000 yrs V s hock = 5000 km s _1 , and diam- 
eter of the remnant after 1000 yrs D = 17 pc. 
The initial magnetic field configuration is that 
suggested by Bocchino et al. (201 1) in the case 
of SN 1006 and resulting from the comparison 
of radio observations with MHD models. The 
initial radius of the remnant is 0.5 pc. We fol- 
low the remnant evolution for 1000 yrs. Special 
emphasis was placed on capturing the enor- 
mous range in spatial scales in the remnant: in 
fact we need to perform simulations with suffi- 
cient spatial resolution to capture the structures 
of the ejecta in the initial system configuration 
and to follow the evolutions of these structures. 
To this end, we exploited the adaptive mesh 
capabilities of the FLASH code by using 11 
nested levels of resolution; the effective mesh 
size was 8192 x 8192 x 8192. 



Orlando et al. 



Ejecta clumping and back-reaction of accelerated cosmic rays in SNRs 





Fig. 1. 3D volume rendering describing the spa- 
tial distribution of the effective adiabatic index for 
a model with an exponential profile of the initial 
ejecta density after 1000 yrs. The index y cff is min- 
imum in red regions (see color bar). The violet sur- 
face tracks the ejecta material and the white lines are 
sampled magnetic field lines. 



3. Results 

We first analyzed the back-reaction of accel- 
erated CRs, by considering models account- 
ing for the shock modification but without ini- 
tial clumping of ejecta. As an example, Fig. 1 
presents the results for a model with an expo- 
nential profile of the initial ejecta density after 
1000 yrs of evolution. In this model we have 
assumed that the back-reaction of accelerated 
CRs is modulated by the obliquity angle, be- 
ing more efficient at parallel shocks where the 
effective adiabatic index is minimum, namely 
Teff = 1-1 (we are assuming extreme energy 
losses to accelerate the CRs). 

We investigated the effect of accelerated 
CRs on the separation between the blast wave 
and the contact discontinuity. To this end, we 
derived the azimuthal profiles of the ratio be- 
tween the blast wave and contact discontinuity 
radii Rbw/Rcd from the models and compared 
them with that derived from the observations 
of SN 1006 (Miceli et al. 2009). In general, we 
found that the modeled profiles are higher than 
those observed and are modulated by the obliq- 
uity angle (see black line in Fig. 3). The ob- 
servations can be reproduced only in limited 
regions where the effects of accelerated CRs 




Fig. 2. 3D rendering of mass density for a model 
with an exponential profile of ejecta after 1000 yrs 
of evolution. The initial clumps have a size of the 
order of 4% of the initial remnant radius and a max- 
imum density contrast of 5. The white contour en- 
closes the ejecta material. 



is higher. This result suggests that the obser- 
vations could be reproduced only if the back- 
reaction of accelerated CRs is extreme (i.e. 
y e ff ~ 1.1) at all obliquity angles. 

As a next step, we investigated the effects 
of ejecta clumping on the evolution and mor- 
phology of the remnant, considering models 
without back-reaction of accelerated CRs and 
accounting for ejecta clumping. As an exam- 
ple, Fig. 2 shows the 3D rendering of mass den- 
sity for a model with an exponential profile of 
ejecta. The figure shows that: 1) the RT mix- 
ing reaches the forward shock front perturb- 
ing the remnant outlines, 2) clumps and fila- 
mentary structures are evident within the rem- 
nant, and 3) clumps of ejecta are close to or 
protruding beyond the main blast wave as ob- 
served in SN 1006 by Rakowski et al. (2011). 
In general, we found that increasing the initial 
size of the clumps or increasing their density 
contrast, both the perturbation of the remnant 
outlines and the occurrence of ejecta protru- 
sions increase. Also in this case, we compared 
the azimuthal profiles of the ratio /?bw/^cd de- 
rived from the models with that observed in 
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Fig. 3. Azimuthal profiles of the ratio between 
the forward shock and the contact discontinuity 
radii for models without ejecta clumping and back- 
reaction of accelerated CRs (Ml, red line), with 
ejecta clumping and without accelerated CRs (M2, 
blue), and with accelerated CRs and without ejecta 
clumping (M3, black). All models assume an expo- 
nential profile of the initial ejecta density; M2 as- 
sumes a minimum y e g = 1.1; M3 assumes a maxi- 
mum density contrast 2.5 and an initial size of ejecta 
4% of the remnant radius. The green line marks the 
profile observed in SN 1006 (Miceli et al. 2009). 



SN 1006. We found that the initial clumping 
of ejecta makes the azimuthal profiles Rbw/Rcd 
fairly uniform and lower than expected for 
models without clumping. In particular, we 
found that: the larger the size of initial clumps 
of ejecta, the lower the value of the ratio; the 
higher the initial density contrast, the lower the 
value of the ratio. 

As an example, in Fig. 3 we compare the 
azimuthal profiles derived for models account- 
ing for only one of the effects considered in this 
paper (either back-reaction of accelerated CRs 
or ejecta clumping), and for a model without 
both physical effects. Our analysis has shown 
that the observed profile can be reproduced by 
models with a maximum density contrast of 
ejecta ranging between 2.5 and 5, and an ini- 
tial size of ejecta clumps of the order of 4% of 
the remnant radius (see, for instance, blue line 
in Fig. 3). 

4. Summary and conclusions 

We investigated the role of ejecta clumping and 
back-reaction of accelerated CRs on the evo- 



lution of young SNRs. To this end, we have 
developed a 3D MHD model including con- 
sistently the back-reaction of accelerated CRs 
and the ejecta clumping. Particular attention 
has been devoted to the spatial resolution that 
is necessary to describe appropriately the struc- 
ture of the ejecta; to this end, we have exploited 
the AMR capabilities of the FLASH code. Our 
starting point was the observational evidence 
that the azimuthal profile of the ratio between 
the forward shock and the contact discontinuity 
radii is fairly uniform and lower than expected 
for a non-modified shock. 

We found that the back-reaction of acceler- 
ated CRs alone cannot reproduce the observa- 
tions unless the CRs energy losses are extreme 
(that is the effective adiabatic index should be 
of the order of 1.1) and independent on the 
obliquity angle. We also found that the clump- 
ing of ejecta is an important factor to reproduce 
the observed values of the ratio and its obliq- 
uity dependence. We conclude that the ejecta 
clumping is important in the description of the 
evolution of young SNRs. 
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